if(exists("snakemake")) {
input_list <- snakemake@input
output_list <- snakemake@output
} else {
# inputs:
input_list <- list(
oregon_rivers = "inputs/Rivers_OR/Rivers_OR.shp"
)
# outputs:
output_list <- list(
texmap = "tex/images/map-crop.pdf"
)
}
# add the other entries based on the final texmap
output_list$map = "results/map_of_samples/map.pdf"
output_list$cropmap = "results/map_of_samples/map-crop.pdf"
# we create the necessary output directories like this:
dump <- lapply(output_list, function(x)
dir.create(dirname(x), recursive = TRUE, showWarnings = FALSE)
)
library(terra)
library(sf)
library(tidyverse)
library(ggspatial)
#library(ggsn)
library(cowplot)
library(plotly)
library(maps) # needed for map_data()
library(mapproj) # making it explicit for renv
if(!file.exists("geo-spatial/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif")) {
dir.create("geo-spatial", showWarnings = FALSE)
download.file(
url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/raster/HYP_HR_SR_OB_DR.zip",
destfile = "geo-spatial/HYP_HR_SR_OB_DR.zip"
)
unzip(
zipfile = "geo-spatial/HYP_HR_SR_OB_DR.zip",
exdir = "geo-spatial/HYP_HR_SR_OB_DR"
)
file.remove("geo-spatial/HYP_HR_SR_OB_DR.zip")
}
if(!file.exists("geo-spatial/ne_10m_admin_1_states_provinces_lines/ne_10m_admin_1_states_provinces_lines.shp")) {
dir.create("geo-spatial", showWarnings = FALSE)
download.file(
url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces_lines.zip",
destfile = "geo-spatial/ne_10m_admin_1_states_provinces_lines.zip"
)
unzip(
zipfile = "geo-spatial/ne_10m_admin_1_states_provinces_lines.zip",
exdir = "geo-spatial/ne_10m_admin_1_states_provinces_lines"
)
file.remove("geo-spatial/ne_10m_admin_1_states_provinces_lines.zip")
}
if(!file.exists("geo-spatial/ne_10m_coastline/ne_10m_coastline.shp")) {
dir.create("geo-spatial", showWarnings = FALSE)
download.file(
url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip",
destfile = "geo-spatial/ne_10m_coastline.zip"
)
unzip(
zipfile = "geo-spatial/ne_10m_coastline.zip",
exdir = "geo-spatial/ne_10m_coastline"
)
file.remove("geo-spatial/ne_10m_coastline.zip")
}
This is from California Nat Res Agency. Find it at: https://data.cnra.ca.gov/dataset/national-hydrography-dataset-nhd
# note. Once when I did this, download.file did not see capable of downloading this file on my Mac
# (the download is super sloooooow), and that time I
# ended up downloading it with Chrome. But it seems those connection problems
# had been fixed when I did this again.
if(!file.exists("geo-spatial/NHD_Major_Rivers_and_Creeks/Major_Rivers_and_Creeks.shp")) {
dir.create("geo-spatial", showWarnings = FALSE)
download.file(
url = "https://data.cnra.ca.gov/dataset/511528b2-f7d3-4d86-8902-cc9befeeeed5/resource/7d1e7e44-81b1-43fe-95f6-1862eea6ac24/download/nhd_major_rivers_and_creeks.zip",
destfile = "geo-spatial/nhd_major_rivers_and_creeks.zip"
)
unzip(
zipfile = "geo-spatial/nhd_major_rivers_and_creeks.zip",
exdir = "geo-spatial"
)
file.remove("geo-spatial/nhd_major_rivers_and_creeks.zip")
}
nat.earth <- rast("geo-spatial/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif") #
state_prov <- st_read("geo-spatial/ne_10m_admin_1_states_provinces_lines/ne_10m_admin_1_states_provinces_lines.shp")
## Reading layer `ne_10m_admin_1_states_provinces_lines' from data source
## `/Users/eriq/Documents/git-repos/california-chinook-microhaps/geo-spatial/ne_10m_admin_1_states_provinces_lines/ne_10m_admin_1_states_provinces_lines.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 10179 features and 57 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -178.1371 ymin: -49.25087 xmax: 178.4486 ymax: 81.12853
## Geodetic CRS: WGS 84
coastline <- st_read("geo-spatial/ne_10m_coastline/ne_10m_coastline.shp")
## Reading layer `ne_10m_coastline' from data source
## `/Users/eriq/Documents/git-repos/california-chinook-microhaps/geo-spatial/ne_10m_coastline/ne_10m_coastline.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4133 features and 3 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -180 ymin: -85.22194 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
all_rivers <- st_read("geo-spatial/NHD_Major_Rivers_and_Creeks/Major_Rivers_and_Creeks.shp") %>%
st_zm() %>%
st_transform(., st_crs(state_prov))
## Reading layer `Major_Rivers_and_Creeks' from data source
## `/Users/eriq/Documents/git-repos/california-chinook-microhaps/geo-spatial/NHD_Major_Rivers_and_Creeks/Major_Rivers_and_Creeks.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 220702 features and 17 fields
## Geometry type: LINESTRING
## Dimension: XYZM
## Bounding box: xmin: -124.4252 ymin: 32.50998 xmax: -113.9205 ymax: 42.6784
## z_range: zmin: -2e-04 zmax: 1476.791
## m_range: mmin: -1.797693e+308 mmax: 100
## Geodetic CRS: NAD83 + NAVD88 height
or_rivers <- st_read(input_list$oregon_rivers)
## Reading layer `Rivers_OR' from data source
## `/Users/eriq/Documents/git-repos/california-chinook-microhaps/inputs/Rivers_OR/Rivers_OR.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 8 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -124.4226 ymin: 38.65888 xmax: -108.1906 ymax: 49.27932
## Geodetic CRS: WGS 84
# important to put them in this order and named like this
domain <- c(
xmin = -127.9,
xmax = -117.4,
ymin = 36,
ymax = 43.25
)
sf_use_s2(FALSE)
nat_crop <- crop(nat.earth, y = ext(domain))
state_subset <- st_crop(state_prov, domain)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
coastline_cropped <- st_crop(coastline, domain)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
We plot this with plotly so we can see the names of the rivers easily, and whittle it down to the ones that we need.
# I don't know the gnis_id for mill creek, and there are a lot of Mill Creeks,
# so lets get it:
mill_area <- c(
xmin = -121.4,
xmax = -122.3,
ymin = 40,
ymax = 40.3
)
mill_candi <- all_rivers %>%
filter(str_detect(gnis_name, "Mill Creek")) %>%
st_crop(mill_area)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# same with Blue Creek
blue_area <- c(
xmin = -124,
xmax = -123.5,
ymin = 41.2,
ymax = 41.8
)
blue_candi <- all_rivers %>%
filter(str_detect(gnis_name, "Blue Creek")) %>%
st_crop(blue_area)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# rogue river:
rogue <- or_rivers %>%
filter(name == "Rogue") %>%
mutate(gnis_name = "Rogue River")
our_rivers <- all_rivers %>%
filter(gnis_name %in% c(
"Sacramento River",
"San Joaquin River",
"Feather River",
"Russian River",
"Trinity River",
"Klamath River",
"Eel River",
"Smith River"
#"Deer Creek",
#"Battle Creek"
) |
gnis_id %in% c(
"00220293", "00237425", "00263498", "00266522", # these are butte creek
"00233775", "01655075" # these are deer
# "00218740", "00229640", "00234966" # these are battle creek
)
) %>%
bind_rows(mill_candi) %>%
bind_rows(blue_candi) %>%
bind_rows(rogue)
# get candidate positions for notations
labels <- read_tsv("inputs/map-notations.tsv", comment = "#")
g <- ggplot() +
geom_sf(data = our_rivers, aes(colour = gnis_name)) +
geom_sf(data = coastline_cropped) +
geom_point(data = labels, aes(x = label_long, y = label_lat), colour = "red") +
geom_point(data = labels, aes(x = tip_long, y = tip_lat), colour = "blue", size = 0.2)
ggplotly(g)
That looks good. So let us proceed.
base_map <- ggplot() +
ggspatial::layer_spatial(nat_crop) +
geom_sf(data = state_subset, color = "gray30", fill = NA) +
geom_sf(data = coastline_cropped, color = "gray30", fill = NA, linewidth = 0.15) +
geom_sf(data = our_rivers, colour = "blue", linejoin = "round", lineend = "round")
Now add stuff to that:
bit <- 0.0
kick <- 0.4
source("R/colors.R")
labels_1line <- labels %>%
group_by(map_text) %>%
filter(rank == max(rank)) %>%
ungroup()
# define the shapes
point_shapes <- c(
`Fall run` = 21,
`Spring run` = 22,
`Winter run` = 23,
`Late-fall run` = 24
)
mm <- base_map +
geom_segment(data = labels, aes(x = tip_long, y = tip_lat, xend = label_long, yend = label_lat), colour = "black", linewidth = 0.4) +
geom_point(data = labels, aes(x = label_long + bit + (1 - hjust * 2) * kick * (rank - 1), y = label_lat, fill = run_timing, shape = run_timing), size = 3.5) +
geom_label(
data = labels_1line,
aes(
x = label_long + bit + (1 - hjust * 2) * kick * ( 0.65 + rank - 1),
y = label_lat,
label = map_text,
hjust = hjust
),
size = 2.2,
label.padding = unit(0.09, "lines"),
lineheight = 0.85
) +
scale_fill_manual(values = run_time_colors) +
scale_shape_manual(values = point_shapes) +
theme_bw() +
theme(
panel.border = element_rect(colour = "black", linewidth = 1),
axis.text.x = element_text(size = 8, family = "serif", angle = 35, hjust = 1),
axis.text.y = element_text(size = 8, family = "serif"),
axis.title.y = element_text(family = "serif", size = 10),
axis.title.x = element_text(family = "serif", vjust = 2, size = 10),
plot.margin = margin(0, 0.1, 0, 0.15, "cm"),
legend.position = "none"
) +
xlab("Longitude") +
ylab("Latitude") +
coord_sf(
expand = FALSE,
) +
guides(fill = guide_legend(title = "Run Timing:", nrow = 2)) +
ggspatial::annotation_north_arrow(location = "br", height = unit(10, "mm"), width = unit(5, "mm"), style = north_arrow_fancy_orienteering()) +
ggspatial::annotation_scale(location = "br", height = unit(0.1, "cm"))
wrld <- map_data("world")
states <- map_data("state")
domain_df <- tibble(point = 1:length(domain), long = rep(domain[1:2], each = 2), lat = c(domain[3:4], rev(domain[3:4])))
inset_world <- ggplot() +
geom_polygon(data = wrld, aes(x = long, y = lat, group = group), colour = "black", fill = "gray90", linewidth = 0.1) +
geom_path(data = states, aes(x = long, y = lat, group = group), colour = "black", linewidth = 0.1) +
geom_polygon(data = domain_df, mapping = aes(x = long, y = lat), colour = "red", fill = "red", alpha = 0.3, linewidth = 0.1) +
coord_map("ortho", orientation = c(30, -114, 0), xlim = c(-150, -75), ylim = c(28, 60)) +
theme_bw() +
labs(x = NULL, y = NULL) +
theme(
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(0, 0, -1, -1), "mm")
)
inset_world
final_map <- ggdraw() +
draw_plot(mm) +
draw_plot(inset_world, x = 0.215, y = 0.139, width = 0.20, height = 0.17)
ggsave(final_map, filename = output_list$map, width = 5, height = 3.5)
# let's crop that down for use, too
CALL <- paste("pdfcrop ", output_list$map, collapse = " ")
system(CALL)
file.copy(from = output_list$cropmap, to = output_list$texmap, overwrite = TRUE)
## [1] TRUE
We can do this and just kick them directly into the tex/images folder.
rc_col_tib <- enframe(run_time_colors) %>%
mutate(
nospace = str_replace_all(name, " ", "-"),
file = str_c("tex/images/", nospace, "-ball.pdf")
)
for(r in 1:nrow(rc_col_tib)) {
tmp <- rc_col_tib[r,]
g <- ggplot(tmp, aes(x = 1, y = 1)) +
geom_point(shape = point_shapes[tmp$name], fill = tmp$value, size = 3.5) +
theme_void()
ggsave(g, filename = tmp$file)
CALL <- str_c("pdfcrop ", tmp$file)
system(CALL)
file.remove(tmp$file)
}
sessioninfo::session_info()
## ─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.4.3 (2025-02-28)
## os macOS Monterey 12.7.5
## system aarch64, darwin20
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Denver
## date 2025-03-27
## pandoc 3.6.4 @ /Users/eriq/Documents/git-repos/california-chinook-microhaps/.snakemake/conda/def120e971fb2c87a8c042b4c2c09bdb_/bin/ (via rmarkdown)
## quarto 1.4.550 @ /usr/local/bin/quarto
##
## ─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## ! package * version date (UTC) lib source
## P bit 4.6.0 2025-03-06 [?] CRAN (R 4.4.1)
## P bit64 4.6.0-1 2025-01-16 [?] CRAN (R 4.4.1)
## P bslib 0.9.0 2025-01-30 [?] CRAN (R 4.4.1)
## P cachem 1.1.0 2024-05-16 [?] CRAN (R 4.4.1)
## P class 7.3-23 2025-01-01 [?] CRAN (R 4.4.3)
## P classInt 0.4-11 2025-01-08 [?] CRAN (R 4.4.1)
## P cli 3.6.4 2025-02-13 [?] CRAN (R 4.4.1)
## P codetools 0.2-20 2024-03-31 [?] CRAN (R 4.4.3)
## P colorspace 2.1-1 2024-07-26 [?] CRAN (R 4.4.1)
## P cowplot * 1.1.3 2024-01-22 [?] CRAN (R 4.4.0)
## P crayon 1.5.3 2024-06-20 [?] CRAN (R 4.4.1)
## P crosstalk 1.2.1 2023-11-23 [?] CRAN (R 4.4.0)
## P data.table 1.17.0 2025-02-22 [?] CRAN (R 4.4.1)
## P DBI 1.2.3 2024-06-02 [?] CRAN (R 4.4.1)
## P digest 0.6.37 2024-08-19 [?] CRAN (R 4.4.1)
## P dplyr * 1.1.4 2023-11-17 [?] CRAN (R 4.4.0)
## P e1071 1.7-16 2024-09-16 [?] CRAN (R 4.4.1)
## P evaluate 1.0.3 2025-01-10 [?] CRAN (R 4.4.1)
## P farver 2.1.2 2024-05-13 [?] CRAN (R 4.4.1)
## P fastmap 1.2.0 2024-05-15 [?] CRAN (R 4.4.1)
## P forcats * 1.0.0 2023-01-29 [?] CRAN (R 4.4.0)
## P generics 0.1.3 2022-07-05 [?] CRAN (R 4.4.1)
## P ggplot2 * 3.5.1 2024-04-23 [?] CRAN (R 4.4.0)
## P ggspatial * 1.1.9 2023-08-17 [?] CRAN (R 4.4.0)
## P glue 1.8.0 2024-09-30 [?] CRAN (R 4.4.1)
## P gtable 0.3.6 2024-10-25 [?] CRAN (R 4.4.1)
## P hms 1.1.3 2023-03-21 [?] CRAN (R 4.4.0)
## P htmltools 0.5.8.1 2024-04-04 [?] CRAN (R 4.4.1)
## P htmlwidgets 1.6.4 2023-12-06 [?] CRAN (R 4.4.0)
## P httr 1.4.7 2023-08-15 [?] CRAN (R 4.4.0)
## P jquerylib 0.1.4 2021-04-26 [?] CRAN (R 4.4.0)
## P jsonlite 1.9.1 2025-03-03 [?] CRAN (R 4.4.1)
## P KernSmooth 2.23-26 2025-01-01 [?] CRAN (R 4.4.3)
## P knitr 1.49 2024-11-08 [?] CRAN (R 4.4.1)
## P labeling 0.4.3 2023-08-29 [?] CRAN (R 4.4.1)
## P lazyeval 0.2.2 2019-03-15 [?] CRAN (R 4.4.1)
## P lifecycle 1.0.4 2023-11-07 [?] CRAN (R 4.4.1)
## P lubridate * 1.9.4 2024-12-08 [?] CRAN (R 4.4.1)
## P magrittr 2.0.3 2022-03-30 [?] CRAN (R 4.4.1)
## P mapproj * 1.2.11 2023-01-12 [?] RSPM
## P maps * 3.4.2.1 2024-11-10 [?] RSPM
## P munsell 0.5.1 2024-04-01 [?] CRAN (R 4.4.1)
## P pillar 1.10.1 2025-01-07 [?] CRAN (R 4.4.1)
## P pkgconfig 2.0.3 2019-09-22 [?] CRAN (R 4.4.1)
## P plotly * 4.10.4 2024-01-13 [?] CRAN (R 4.4.0)
## P proxy 0.4-27 2022-06-09 [?] CRAN (R 4.4.1)
## P purrr * 1.0.4 2025-02-05 [?] CRAN (R 4.4.1)
## P R6 2.6.1 2025-02-15 [?] CRAN (R 4.4.1)
## P ragg 1.3.3 2024-09-11 [?] CRAN (R 4.4.1)
## P Rcpp 1.0.14 2025-01-12 [?] CRAN (R 4.4.1)
## P readr * 2.1.5 2024-01-10 [?] CRAN (R 4.4.0)
## renv 1.1.4 2025-03-20 [1] CRAN (R 4.4.1)
## P rlang 1.1.5 2025-01-17 [?] CRAN (R 4.4.1)
## P rmarkdown 2.29 2024-11-04 [?] CRAN (R 4.4.1)
## P sass 0.4.9 2024-03-15 [?] CRAN (R 4.4.0)
## P scales 1.3.0 2023-11-28 [?] CRAN (R 4.4.0)
## P sessioninfo 1.2.3 2025-02-05 [?] CRAN (R 4.4.1)
## P sf * 1.0-20 2025-03-24 [?] RSPM (R 4.4.0)
## P stringi 1.8.4 2024-05-06 [?] CRAN (R 4.4.1)
## P stringr * 1.5.1 2023-11-14 [?] CRAN (R 4.4.0)
## P systemfonts 1.2.1 2025-01-20 [?] CRAN (R 4.4.1)
## P terra * 1.8-29 2025-02-26 [?] CRAN (R 4.4.1)
## P textshaping 1.0.0 2025-01-20 [?] CRAN (R 4.4.1)
## P tibble * 3.2.1 2023-03-20 [?] CRAN (R 4.4.0)
## P tidyr * 1.3.1 2024-01-24 [?] CRAN (R 4.4.1)
## P tidyselect 1.2.1 2024-03-11 [?] CRAN (R 4.4.0)
## P tidyverse * 2.0.0 2023-02-22 [?] CRAN (R 4.4.0)
## P timechange 0.3.0 2024-01-18 [?] CRAN (R 4.4.1)
## P tzdb 0.4.0 2023-05-12 [?] CRAN (R 4.4.0)
## P units 0.8-7 2025-03-11 [?] CRAN (R 4.4.1)
## P vctrs 0.6.5 2023-12-01 [?] CRAN (R 4.4.0)
## P viridisLite 0.4.2 2023-05-02 [?] CRAN (R 4.4.1)
## P vroom 1.6.5 2023-12-05 [?] CRAN (R 4.4.0)
## P withr 3.0.2 2024-10-28 [?] CRAN (R 4.4.1)
## P xfun 0.51 2025-02-19 [?] CRAN (R 4.4.1)
## P yaml 2.3.10 2024-07-26 [?] CRAN (R 4.4.1)
##
## [1] /Users/eriq/Documents/git-repos/california-chinook-microhaps/renv/library/macos/R-4.4/aarch64-apple-darwin20
## [2] /Users/eriq/Library/Caches/org.R-project.R/R/renv/sandbox/macos/R-4.4/aarch64-apple-darwin20/f7156815
##
## * ── Packages attached to the search path.
## P ── Loaded and on-disk path mismatch.
##
## ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Running the code and rendering this notebook required approximately this much time:
Sys.time() - start_time
## Time difference of 16.5015 secs